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A modified Kuramoto model of synchronization in a finite discrete system of 
locally coupled oscillators is studied. The model consists of N oscillators with 
random natural frequencies arranged on a ring. It is shown analytically and nu- 
merically that finite-size systems may have many different synchronized stable 
solutions which are characterised by different values of the winding number. The 
lower bound for the critical coupling kc is given, as well as an algorithm for its 
exact calculation. It is shown that in general phase-locking does not lead to phase 
coherence in ID. 

PACS numbers: 05.45.Xt 

1. Introduction 

The phenomenon of collective synchronization, in which a large population of 
elements oscillating with different frequencies spontaneously locks to a common 
frequency, is the subject of extensive research in physics, biology, chemistry, and 
social sciences. Biological examples at various levels of complexity contain pace- 
maker cells in the heart [1,2], the metabolic synchrony in yeast cell suspensions 
Em or the synchronously flashing swarm of fireflies [5, 6J. The studies of syn- 
chronizing man-made systems include arrays of lasers [7, 8l or superconducting 
Josephson junctions 1 9,. 10]. 
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The Kuramoto model fTP|, which treats the system as an ensemble of limit- 
cycle oscillators described only by their phases, proved to be a very successful 
approach to the problem of synchronization. The natural frequencies of the ele- 
ments of the population are drawn from some prescribed distribution. It has been 
shown analytically that the system with a mean-field coupling exhibits a phase 
transition in the thermodynamic limit: if the coupling exceeds a certain threshold, 
a macroscopic set of the oscillators spontaneously synchronize. 

On the other hand, as shown by Strogatz and MiroUo [12], an infinite ID sys- 
tem does not synchronize. However, in real-life application only finite systems 
appear. It is, therefore, interesting to see whether finite ID Kuramoto systems may 
synchronize. 



2. A Local One-Dimensional Kuramoto Model 

The topology of interactions of the studied system forms a ring, i.e. we as- 
sume the lattice to be one-dimensional with nearest-neighbour coupling, and peri- 
odic boundary conditions. The local Kuramoto model of synchronizing limit-cycle 
oscillators is given by the equations of motion (first-order, nonlinear ODEs) 

ei(0 = (Oi-^[sin(e,-(0-e/-i(0) + sin(ei(0-e,+i(0)] , (i) 

for / = 1, . . . , A'^, where 6, are the phases of the oscillators, co, are the natural fre- 
quencies taken from a certain distribution, and k is the coupling constant. 

3. The form of stable synchronized solutions 

We look for stationary synchronized solutions given by 

Qi{t) = ^t + (Sfi ,i = \,...,N , (2) 

where Q. is the common frequency, and the set of N constant phases {(|)i}i=i,...,Ar 
characterises a given solution. The insertion of this condition into the equations of 
motion ([T|l yields a linear system of equations for sines of constant phase differ- 
ences: 

2 2 

As= -(0}-ae) = -A , (3) 
k k 

where e= [1, 1, . . . , 1]^, s= [sin((t)i -(t)A,),sin((t)2 -(t)i), . . . ,sin((|)Ar -(t)A?_i)], (0 = 
[coi,a)2, .-,(Ow]^, A= [5i,52,...,5a?]^ = [coi -Q.,0}2-^, --,0)n -^V are the 
deviations of natural frequencies from the mean frequency, and the matrix A is 
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^iven by 



(4) 



1 



Since the vector e = [1 , 1 , . . . , 1]^ is an eigenvector of A to the zero eigenvalue, 
the determinant of the matrix vanishes. Consequently, the linear equation ([3]l has 
one free parameter p € [—1,1] (it behaves as en element of the vector s, hence the 
constraint), which allows the solutions to appear, as shown below. 

The system is easily solvable but leads to some constraints on what systems 
(e.g. having certain distributions of natural frequencies) can synchronize, which 
comes from the observation that \si\ < 1. The condition of solvability of the linear 
equation leads to the observation that the synchronized frequency is equal to the 
mean of the distribution of natural frequencies (the mean always exists in a finite 
population) 

1 ^ 

a = - y (0, . (5) 

The synchronized solutions take the form of a set of phase differences between 
the neighbouring oscillators 



N-i 



- ^i-i = arcsin P + t XI , / = 1, . . . , 



(6) 



where G [— 1, 1] is a parameter of the solution of the linear equation ([3]l. For 
i = N we assume that the sum equals 0, i.e. — <\)n-\ = arcsin p. 

Above, it is assumed that all inverse functions of sines Si = sin (([); — 
should be taken as (|),- — (j),_i = arcsin 5,- e (— 7l/2,7l/2] for the sake of stability of 
the solution (as discussed in Sec. 8 in more detail). For any of — 1 independent 
phase differences it is possible to take another inverse, n — arcsin^,- e (7i/2,37i/2], 
therefore in total one can obtain 2^^^ possibilities of how the solution can look 
like. 



4. The number of stable solutions 

For a given coupling k and frequency distribution one finds stable solutions 
by solving for p £ [—1, 1] the equation obtained from summing up the N phase 
differences (|6]) 



N-l 



2 



X arcsin + t 51 '^y ~ ^'^'^ ' (7) 
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Fig. 1. The illustration of the equation (jT) for = 6. Intersection of the curves (left-hand 
side of the equation) with the 0,±27t lines means that 0, 1, 2 or 3 solutions appear. The 
broken lines indicate limiting values obtained by the left-hand side. 




Fig. 2. The idea of the winding number: the circle represents phase differences (|), — £ 
[0,27t). The digits are oscillators' indices. The right-hand solution has m = 1, all the phase 
differences are positive and sum up to 2%. The left-hand has m = 0, one of the differences 
is negative and the sum is zero. 



where m = - [A^/4J , - [A^/4J + 1 , . . . , - 1 , 0, 1 , . . . , [A^/4J , which we call the "wind- 
ing number" ([.J denotes "'floor"', i.e. the largest integer number which is lower 
than or equal to the given number). Thus, there can exist more than one stable so- 
lution. Indeed, there can be at most 1 + 2 • [N/4\ synchronized solutions differing 
in winding numbers. 
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5. Restrictions on synchronizing systems 

The form of interactions results in the fact that for synchronization to appear, 
all sums of deviations from the mean frequency must obey the inequality 



<1 , (8) 



which binds the coupling to the natural frequency distribution. In the limit N oo 
independently of the frequency distribution, as shown by Strogatz and Mirollo 
lfT2l . the behaviour of the sums is given by 

ri/2 



max I V A,| ~A^i/2 _ 

l<i<N f-'^ 

It comes from the fact that walking along the ring we can think of the subsequent 
random A, as the steps of a one-dimensional random walk. This result was first 
obtained for Gaussian distribution but holds for any distribution of independently 
distributed natural frequencies. 

It follows that if the critical coupling constant kc obeys the inequality (O, it is 
infinite in the limit N —foo. Consequently, infinite systems cannot synchronize, and 
so we restrict ourselves to the study of finite systems. The conclusion is compatible 
with the calculation of lower critical dimension by Daido iflBl . 

The behaviour resembles that of the Ising model, which does not exhibit phase 
transition in one dimension, either (a brief discussion on the Ising model's solu- 
tion on a ring as well as a commentary on the lack of phase transition in ID can 
be found in fT4l). In the Ising model the phase transition occurs neither in infinite 
nor finite systems, while the one-dimensional Kuramoto model exhibits synchro- 
nization for a finite number of oscillators. 



6. Phase vs. frequency synchronization 

The order parameter for the original Kuramoto model, 

rit)e''^^'^ = l:t'"'^'^ ' (10) 

indicates how well synchronized the oscillators are: if r = 0, the system is totally 
incoherent, and if r = 1 , it is fully synchronized. It defines the critical couphng 
constant fc^ for which the continuous phase transition occurs: r = only for k <kc, 
otherwise r > 0. 

In the local model the given definition does not serve its purpose. r{k) pro- 
vides the information on phase coherence which is a stronger condition than just 
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Fig. 3. The order parameter r{k). The synchronized solutions can be phase incoherent, 
hence r « for m = ± 1 . 



frequency synchronization. Only the solution having a zero winding number can 
take values r > in the limit of large k. It results from the phase differences taking 
uniform values when k increases (which can be seen already in the equation Q). 



7. The critical coupling 

It is interesting to find the critical coupling, where the phase transitions takes 
place for a given distribution of frequencies. Here, as we restrict ourselves to 
finite-size systems, just indicates the strength of interactions for which the first 
synchronized solution appears. To obtain kc for a given system (i.e. size N and 
frequencies) one needs to calculate min = | min,=o,...,A'-i (L;=i'^;) + 1> max = 
|max,=o N-i (L;=i'^;) ~ 1 arid next determine the values assumed by the left- 
hand side of the equation Q 



^ arcsin ij^^j — min j , ^ arcsin ij^^j — max 



(11) 



The smallest value of k for which one of the ends of the above interval equals Imn 
is the kc. This is tantamount to finding when (in terms of k) solutions appear, and 
choosing the first one. 



8. The stability of solutions 

We analyse the linear perturbation of the stationary solutions. Let {6;(0}£i = 
{Q.t + <^i}f^i be a synchronized solution of the system ([T]). We perturb it V/ : 
9,(f) Qi{t) +Ui{t), where \ui{t)\ <C 1, and substitute it into the equations of 
motion ([B. Having linearised the right-hand side with respect to m, we obtain the 



J.Ochab,P.F.Gora printed on August 31, 2009 



7 



evolution equation of the perturbation 



C\+C2 -C2 
-C2 C2 + C3 -C3 
-C3 C3 + C4 



-C4 



-CN CN + Cl 



where m = [ui{t),U2it), ■ . .,UN{t)Y,Ci = cos ((]),•- (t);_i). 

We consider the phase differences from the interval (j),- — (|);_ i G [—71, 7i] (mod 2%) . 
Due to its cyclic structure, the matrix in ([T2l ) always has a zero eigenvalue, so 
that every solution is neutrally stable with respect to the homogeneous translations 
("rigid rotations" of the system): V/ : ([),• (]); + All the other eigenvalues depend 
on the values of c,. 

We use the Gershgorin theorem (which can be found in |[T5l[T6l ). which states 
that for a matrix ^ all the eigenvalues lie in at least one of the circles 

(so-called Gershgorin discs) given by {z : |z — an \ < /?,}, where 7?, = Y!j=\j^i I^'O' !> 
i=\,...,N. 

Using the theorem to localise the eigenvalues of the perturbation matrix we 
classify the solutions according to their phase differences as follows: 

(i) V/ = 1, . . . ,A'^ : !([); — 1 < 7i/2 => stable solutions (the eigenvalues Xi obey 

ReX; < 0) 

(ii) V/ = 1, . . . ,A'^ : — (j);-! I > 7i/2 (mod27i) => unstable solutions (the eigen- 

values Xi obey ReX,- > 0) 

(iii) 3/ : [(]),• — | < 7i/2 A 3j : — | > '!i/2=> we cannot determine the 

stability of the solution 

The third condition describes the situation where at least one of Gershgorin 
discs intersects both the right and the left imaginary half-planes. The theorem 
does not tell us in which part of the disc the eigenvalues are localised, and whether 
any eigenvalue lies in it at all. In other words, there can exist stable solutions 
outside of (/) and there can exist unstable solutions outside of (//). 



9. The size of basins of attraction 

The numerical evidence indicates that the stable solutions depend on initial 
conditions. The basin of attraction is a connected set containing the given solution. 
If it were otherwise, i.e. if the system reached a given stationary solution starting 
from a set of initial conditions which did not contain the solution, its trajectory 
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Fig. 4. The basins of attraction and the eigenvalues of the perturbation matrix for solutions 
differing in winding number Around k = 0.3 the m = 1 solution appears, then the one 
with m~Q takes over most of its basin of attraction, and at last, near ^ = 2.5 the solution 
m — —\ appears building up its own basin with volume comparable to m— 1 . 



would must have intersected a basin of attraction of another solution. As every 
point of the trajectory is identical to some initial condition, the system should have 
reached the other attractor, which is a contradiction. 

We assume correspondence between the volume of the basins of attraction 
and the stability of solutions, that is the solutions having the maximal negative 
eigenvalue of the perturbation matrix closer to zero (i.e. they are "less stable") 
should have smaller basin of attraction than the "more stable" solutions whose all 
eigenvalues are strongly negative. It arises from the fact that for the system of first- 
order ODEs setting any initial condition non-identical to a given stationary solution 
is in fact perturbation of the solution. Since the negative eigenvalues indicate the 
speed of return to the solution, they also tell us which initial condition leads to 
which attractor. Of course, in the stability analysis we assume the perturbation to 
be small \ui{t) \ ^ 1, so the given argument should be treated as an assumption. 

The simulations suggest that the value of winding number strongly correlates 
with both the eigenvalues and the volume of basins of attraction. The only analyt- 
ical result states that the solutions with large m are marginalised, because of the 
eigenvalues' behaviour 

Xi ~ — ^cos ^-^m^ , (13) 

where m = - [A^/4J , - [A^/4J + 1, . . . , -1,0, 1, . . . , LA^/4J. It is derived from the 
assumption that all the phase differences are equal, which is approximately the 
case for large m and in the limit of large k for any m. 

10. Conclusions 

As has been shown here, the one-dimensional model exhibits synchronization 
only for finite sizes of the system. The stationary synchronized solutions can be 



J.Ochab,P.F.Gora printed on August 3 1 , 2009 



9 



found explicitly, they depend on initial conditions, and their number is proportional 
to A^. The size of the basins of attraction depends on the stability of the attracting 
solutions and, as we have shown, there is a coimection between the winding num- 
ber m and the volume of the basin of attraction, according to which solutions with 
high m are marginalised. 

The order parameter for the Kuramoto model is not very useful here, as it 
indicates only phase coherence, while frequency synchronization may occur here 
without it. The critical coupling constant kc (which in finite systems considered 
here is understood as the smallest coupling for which any synchronized solution 
appears) can easily be computed. 

We hope that the analytically investigated behaviour of the one-dimensional 
model may give a clue, what to expect in the higher-dimensional local models of 
synchronizing oscillators. 
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